This tutorial follows from our introduction to multivariate linear models, extending it by using multivariate linear mixed models.
Useful packages:
library(MCMCglmm)
library(lme4)
library(brms)
library(tidyr)
library(dplyr)
library(corrplot)
library(broom.mixed)
library(dotwhisker)
library(ggplot2); theme_set(theme_bw())
Get data:
data_url <- "http://datadryad.org/bitstream/handle/10255/dryad.8377/dll.csv"
if (!file.exists("dll.csv")) {
download.file(data_url,dest="dll.csv")
}
dll_data <- read.csv("dll.csv")
## make temp a factor (25 vs 30 degrees)
dll_data$temp <- factor(dll_data$temp)
## scale relevant variables (fancier but less repetition than previously)
morph_vars <- c("femur","tibia","tarsus","SCT")
morph_vars_sc <- paste(morph_vars,"s",sep="_")
dll_data2 <- dll_data
## c() drops unwanted structure from the results of scale()
for (i in 1:length(morph_vars)) {
dll_data2[[morph_vars_sc[i]]] <- c(scale(dll_data[[morph_vars[i]]]))
}
The previous expression for the multivariate model was
\[ \mathbf{Y} = \mathbf{XB} + \mathbf{E} \]
In contrast, the expression for the mixed model is
\[ y = \mathbf{X\beta} + \mathbf{Zb} + \epsilon \]
where \(\mathbf{b}\) is a set of Gaussian variables with a variance-covariance matrix \(\mathbf{\Sigma}\) which we estimate.
Suppose we have observations of \(m\) individuals, with \(p\) different observations (traits, or time points, or types of measurement, or …) for each individual. The way we’re going to make this work is to expand (“melt”, or gather() if we’re using tidyr) the data set to be \(mp\) observations long, then treat each individual (which was previously a single row of the data set but now comprises \(p\) rows) as a group (we’ll call this units):
lme4 does not (currently) have a natural syntax for multivariate responses, however, as I eluded to in class, there is an important relationship between multivariate response models and so called “repeated” measures (or longitudinal) models. As such we can use a few tricks to estimate the model in lme4. Below this, I will go through the same model using MCMCglmm, which is a library which has a more natural syntax for such multivariate responses, but is explictly Bayesian, so you need to provide prior distributions.
What we need to first do (for lme4) is to generate a single column that represents the numeric values for our response traits, and then a second variable that stores the trait type.
dll_melt <- (dll_data2
%>% select(-c(femur,tibia,tarsus,SCT)) ## drop unscaled vars
%>% mutate(units=factor(1:n()))
%>% gather(trait,value, -c(units,replicate,line,genotype,temp))
%>% drop_na()
%>% arrange(units)
)
saveRDS(dll_melt, file="dll_melt.rds")
And we can take a look at how this has changed the structure of the data from a “wide” format to a long format
head(dll_data) # original wide
## replicate line genotype temp femur tibia tarsus SCT
## 1 1 line-1 Dll 25 0.590 0.499 0.219 9
## 2 1 line-1 Dll 25 0.550 0.501 0.214 13
## 3 1 line-1 Dll 25 0.588 0.488 0.211 11
## 4 1 line-1 Dll 25 0.588 0.515 0.211 NA
## 5 1 line-1 Dll 25 0.596 0.502 0.207 12
## 6 1 line-1 Dll 25 0.577 0.499 0.207 14
To the long format
head(dll_melt)
## replicate line genotype temp units trait value
## 1 1 line-1 Dll 25 1 femur_s 1.514
## 2 1 line-1 Dll 25 1 tibia_s 0.597
## 3 1 line-1 Dll 25 1 tarsus_s 1.751
## 4 1 line-1 Dll 25 1 SCT_s -1.222
## 5 1 line-1 Dll 25 2 femur_s 0.138
## 6 1 line-1 Dll 25 2 tibia_s 0.654
We can now go ahead and fit the model where we need to include trait as a predictor variable (where each trait is now a “repeated measure” from a particular subject/individual)
t1 <- system.time(
lmer1 <- lmer(value ~ trait:(genotype*temp) - 1 +
(trait-1|line) + (trait-1|units),
data=dll_melt,
control=lmerControl(optCtrl=list(ftol_abs=1e-10),
optimizer="bobyqa",
check.nobs.vs.nlev="ignore",
check.nobs.vs.nRE="ignore"))
)
## Warning in (function (npt = min(n + 2L, 2L * n), rhobeg = NA, rhoend =
## NA, : unused control arguments ignored
-1 to drop intercepttrait:(1+genotype+temp+genotype:temp)trait + trait:genotype + trait:temp + trait:genotype:temp(trait-1|line) means “variance/covariances of traits among lines”-1 so we consider traits, not differences among traits(trait-1|units) specifies “var/cov of traits among units (individuals)”lmer always includes a residual variance term. This is redundant because we have only one data point per individual per trait: lmerControl(...) tells lmer not complainprint(lmer1), summary(lmer1) useful but awkward (16 fixed effect coefficients, we’ll look at graphical summaries insteadfixef() (fixed effects), ranef() (line/indiv deviations from pop mean), coef() (line/indiv-level estimates), VarCorr (random effects var/cov)getME(fit,"theta") extracts themall(abs(getME(lmer1,"theta"))>1e-4) ## check for singularity (OK in this case)
## [1] TRUE
VarCorr(lmer1)
## Groups Name Std.Dev. Corr
## units traitfemur_s 0.654
## traitSCT_s 0.728 0.05
## traittarsus_s 0.565 0.58 0.20
## traittibia_s 0.670 0.91 0.12 0.47
## line traitfemur_s 0.489
## traitSCT_s 0.392 0.14
## traittarsus_s 0.462 0.81 0.37
## traittibia_s 0.473 0.87 0.22 0.83
## Residual 0.462
par(mfrow=c(1,2))
vv1 <- VarCorr(lmer1)
## fix unit variance-covariance by adding residual variance:
diag(vv1$units) <- diag(vv1$units)+sigma(lmer1)^2
corrplot.mixed(cov2cor(vv1$line),upper="ellipse")
corrplot.mixed(cov2cor(vv1$units),upper="ellipse")
Correlations of traits across lines are stronger than correlations within individuals. In both cases correlations are all positive (i.e. first axis of variation is size variation?)
cc1 <- tidy(lmer1,effect="fixed") %>%
tidyr::separate(term,into=c("trait","fixeff"),extra="merge",
remove=FALSE)
dwplot(cc1)+facet_wrap(~fixeff,scale="free",ncol=2)+
geom_vline(xintercept=0,lty=2)
These results tell approximately the same story (coefficients are consistent across traits within a fixed effect, e.g. effect of higher temperature is to reduce scores on all traits).
This is very slow, you probably don’t want to evaluate it on your computer
cc2 <- tidy(lmer1,
effect = "ran_pars",
conf.int = TRUE,
conf.method = "profile",
parallel="multicore",
ncpus=2)
Another option is the MCMCglmm package, which has a natural interface for general multivariate mixed models. It takes a while to get used to the interface, but here is an example. Please check out here for more information about the package, and here for an overview of how to use it.
First I find it easier (given the interface of MCMCglmm) to create a formula with the response variables and predictors. This is only for the fixed effects part of the model.
fmla.MMLM1 <- cbind(femur_s, tibia_s, tarsus_s, SCT_s) ~
trait:(genotype*temp) - 1
Now we need to let MCMCglmm know which family (i.e. distribution) the response variables are. Since all are normal (Gaussian), we can specify it the following way.
fam.test <- rep("gaussian", 4 )
Since MCMCglmm is fundamentally a Bayesian approach, it needs a prior. If you provide no prior by default, it tries a “flat” prior, although this rarely works. In this case I am providing a not-quite-flat prior, but just for the random effects of line and for the residual variances (could also provide them for the fixed effects). Choosing priors for variances and covariances can sometimes be tricky, and for the moment is too much to chew on in this class, so we will use the default prior distributions (inverse-Wishart is what it is called if you really care). If you want to get deep into this, come chat with me. I also have some simple tutorials where I draw out the effects of various common priors for variances and covariances.
prior.model.1 <- list( R = list(V=diag(4)/4, nu=0.004),
G = list(G1=list(V=diag(4)/4, nu=0.004)))
Let’s take a quick look at the prior
prior.model.1
## $R
## $R$V
## [,1] [,2] [,3] [,4]
## [1,] 0.25 0.00 0.00 0.00
## [2,] 0.00 0.25 0.00 0.00
## [3,] 0.00 0.00 0.25 0.00
## [4,] 0.00 0.00 0.00 0.25
##
## $R$nu
## [1] 0.004
##
##
## $G
## $G$G1
## $G$G1$V
## [,1] [,2] [,3] [,4]
## [1,] 0.25 0.00 0.00 0.00
## [2,] 0.00 0.25 0.00 0.00
## [3,] 0.00 0.00 0.25 0.00
## [4,] 0.00 0.00 0.00 0.25
##
## $G$G1$nu
## [1] 0.004
The R matrix prior is for the residuals. This is \(\mathbf{R}_{4,4}\) as we have 4 traits in the VCV matrix. We have have non-zero variances for each trait as the mode for the prior, and 0 for covariances. However this is a weak prior, so even small amounts of data will largely overcome the pull of the prior. Another sensible approach we have used is to make the prior proportional to the overall observed VCV, as there is a large literature on the proportionality of covariance matrices for morphology. This may not be sensible for other types of multivariate response measures.
##,depends.on=c("prior","mod_fm","get_data")}
t2 <- system.time(
MMLM1.fit <- MCMCglmm(fmla.MMLM1,
random=~ us(trait):line,
rcov=~ us(trait):units,
prior= prior.model.1,
data= dll_data2,
family = fam.test,
nitt= 10000, burnin= 2000, thin=10)
)
## Warning: 'cBind' is deprecated.
## Since R version 3.2.0, base's cbind() should work fine with S4 objects
fmla.MMLM1 is the formula object we created abovelmer formulatrait term is a reserved word in MCMCglmm, letting it know we want to fit a multivariate mixed modelMCMCglmm automatically melts the data for us (and assigns the name trait the same way we did manually above)random=~us(trait):line asks MCMCglmm to fit an unstructured covariance matrix for the line term (i.e the different wild type genotypes we are examining).- “Unstructured” means we are estimating the complete 4 x 4 matrix of covariances (= 4*5/2 = 10 elements total)(trait-1|line) in lmer modelMCMCglmm also automatically adds a unitsrcov=~us(trait):units = (trait-1|units)MCMCglmm offers a few other options for variance-covariance structuresThe nitt is how many iterations for the MCMC we want to perform, and the burnin is how many should be ignored at the beginning of the random walk.
The fit takes 22 seconds.
Normally we would spend a fair bit of time on diagnostics of the MCMC, but for now we will just quickly check the trace plots and autocorrelation.
In the fitted object $Sol is for solution, which is the term used for fixed effects in MCMCglmm. $VCV is for the variance-covariance matrix.
library(lattice)
xyplot(MMLM1.fit$Sol[,1:4])
xyplot(MMLM1.fit$Sol[,13:16])
acf(MMLM1.fit$Sol[,1:2])
xyplot(MMLM1.fit$VCV[,1:4])
acf(MMLM1.fit$VCV[,1:3])
Nothing terribly worrying.
summary(MMLM1.fit)
##
## Iterations = 2001:9991
## Thinning interval = 10
## Sample size = 800
##
## DIC: 17451
##
## G-structure: ~us(trait):line
##
## post.mean l-95% CI u-95% CI eff.samp
## traitfemur_s:traitfemur_s.line 0.2965 0.12940 0.488 800
## traittibia_s:traitfemur_s.line 0.2494 0.11156 0.440 800
## traittarsus_s:traitfemur_s.line 0.2268 0.09558 0.396 800
## traitSCT_s:traitfemur_s.line 0.0297 -0.06855 0.148 800
## traitfemur_s:traittibia_s.line 0.2494 0.11156 0.440 800
## traittibia_s:traittibia_s.line 0.2747 0.13994 0.482 800
## traittarsus_s:traittibia_s.line 0.2235 0.10210 0.379 800
## traitSCT_s:traittibia_s.line 0.0468 -0.04449 0.158 800
## traitfemur_s:traittarsus_s.line 0.2268 0.09558 0.396 800
## traittibia_s:traittarsus_s.line 0.2235 0.10210 0.379 800
## traittarsus_s:traittarsus_s.line 0.2651 0.11982 0.439 800
## traitSCT_s:traittarsus_s.line 0.0805 -0.00751 0.203 800
## traitfemur_s:traitSCT_s.line 0.0297 -0.06855 0.148 800
## traittibia_s:traitSCT_s.line 0.0468 -0.04449 0.158 800
## traittarsus_s:traitSCT_s.line 0.0805 -0.00751 0.203 800
## traitSCT_s:traitSCT_s.line 0.1919 0.09279 0.312 585
##
## R-structure: ~us(trait):units
##
## post.mean l-95% CI u-95% CI eff.samp
## traitfemur_s:traitfemur_s.units 0.6441 0.60688 0.6817 800
## traittibia_s:traitfemur_s.units 0.4014 0.36701 0.4355 800
## traittarsus_s:traitfemur_s.units 0.2142 0.18393 0.2413 1100
## traitSCT_s:traitfemur_s.units 0.0257 -0.00674 0.0617 800
## traitfemur_s:traittibia_s.units 0.4014 0.36701 0.4355 800
## traittibia_s:traittibia_s.units 0.6653 0.62308 0.7078 800
## traittarsus_s:traittibia_s.units 0.1799 0.15162 0.2095 897
## traitSCT_s:traittibia_s.units 0.0600 0.02338 0.0918 904
## traitfemur_s:traittarsus_s.units 0.2142 0.18393 0.2413 1100
## traittibia_s:traittarsus_s.units 0.1799 0.15162 0.2095 897
## traittarsus_s:traittarsus_s.units 0.5353 0.50502 0.5723 800
## traitSCT_s:traittarsus_s.units 0.0847 0.05381 0.1153 576
## traitfemur_s:traitSCT_s.units 0.0257 -0.00674 0.0617 800
## traittibia_s:traitSCT_s.units 0.0600 0.02338 0.0918 904
## traittarsus_s:traitSCT_s.units 0.0847 0.05381 0.1153 576
## traitSCT_s:traitSCT_s.units 0.7482 0.70193 0.7935 800
##
## Location effects: cbind(femur_s, tibia_s, tarsus_s, SCT_s) ~ trait:(genotype * temp) - 1
##
## post.mean l-95% CI u-95% CI eff.samp
## traitfemur_s:genotypeDll 0.5062 0.2907 0.7455 699
## traittibia_s:genotypeDll 0.6377 0.4379 0.8731 800
## traittarsus_s:genotypeDll 0.6142 0.4078 0.8096 800
## traitSCT_s:genotypeDll 0.6762 0.4884 0.8625 800
## traitfemur_s:genotypewt 0.1779 -0.0306 0.4093 717
## traittibia_s:genotypewt 0.0170 -0.1994 0.2219 800
## traittarsus_s:genotypewt 0.4418 0.2339 0.6413 800
## traitSCT_s:genotypewt -0.0887 -0.2550 0.1191 800
## traitfemur_s:temp30 -0.8377 -0.9266 -0.7294 800
## traittibia_s:temp30 -0.8091 -0.9134 -0.6873 800
## traittarsus_s:temp30 -1.2379 -1.3353 -1.1509 800
## traitSCT_s:temp30 -0.8386 -0.9588 -0.7244 800
## traitfemur_s:genotypewt:temp30 0.3623 0.2356 0.5186 800
## traittibia_s:genotypewt:temp30 0.4745 0.3148 0.6177 800
## traittarsus_s:genotypewt:temp30 0.5158 0.3968 0.6542 612
## traitSCT_s:genotypewt:temp30 0.7275 0.5759 0.8715 800
## pMCMC
## traitfemur_s:genotypeDll <0.001 **
## traittibia_s:genotypeDll <0.001 **
## traittarsus_s:genotypeDll <0.001 **
## traitSCT_s:genotypeDll <0.001 **
## traitfemur_s:genotypewt 0.093 .
## traittibia_s:genotypewt 0.875
## traittarsus_s:genotypewt <0.001 **
## traitSCT_s:genotypewt 0.333
## traitfemur_s:temp30 <0.001 **
## traittibia_s:temp30 <0.001 **
## traittarsus_s:temp30 <0.001 **
## traitSCT_s:temp30 <0.001 **
## traitfemur_s:genotypewt:temp30 <0.001 **
## traittibia_s:genotypewt:temp30 <0.001 **
## traittarsus_s:genotypewt:temp30 <0.001 **
## traitSCT_s:genotypewt:temp30 <0.001 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Sometimes it is easier to look at the fixed and random effects separately.
summary(MMLM1.fit$Sol)
##
## Iterations = 2001:9991
## Thinning interval = 10
## Number of chains = 1
## Sample size per chain = 800
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## traitfemur_s:genotypeDll 0.5062 0.1132 0.00400 0.00428
## traittibia_s:genotypeDll 0.6377 0.1123 0.00397 0.00397
## traittarsus_s:genotypeDll 0.6142 0.1061 0.00375 0.00375
## traitSCT_s:genotypeDll 0.6762 0.0975 0.00345 0.00345
## traitfemur_s:genotypewt 0.1779 0.1136 0.00402 0.00424
## traittibia_s:genotypewt 0.0170 0.1089 0.00385 0.00385
## traittarsus_s:genotypewt 0.4418 0.1064 0.00376 0.00376
## traitSCT_s:genotypewt -0.0887 0.0974 0.00344 0.00344
## traitfemur_s:temp30 -0.8377 0.0517 0.00183 0.00183
## traittibia_s:temp30 -0.8091 0.0581 0.00206 0.00206
## traittarsus_s:temp30 -1.2379 0.0477 0.00168 0.00168
## traitSCT_s:temp30 -0.8386 0.0598 0.00211 0.00211
## traitfemur_s:genotypewt:temp30 0.3623 0.0720 0.00255 0.00255
## traittibia_s:genotypewt:temp30 0.4745 0.0782 0.00276 0.00276
## traittarsus_s:genotypewt:temp30 0.5158 0.0651 0.00230 0.00263
## traitSCT_s:genotypewt:temp30 0.7275 0.0781 0.00276 0.00276
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## traitfemur_s:genotypeDll 0.278 0.4301 0.5069 0.5845 0.736
## traittibia_s:genotypeDll 0.409 0.5643 0.6364 0.7126 0.857
## traittarsus_s:genotypeDll 0.411 0.5427 0.6149 0.6855 0.813
## traitSCT_s:genotypeDll 0.486 0.6120 0.6781 0.7381 0.858
## traitfemur_s:genotypewt -0.037 0.0939 0.1784 0.2584 0.401
## traittibia_s:genotypewt -0.193 -0.0565 0.0154 0.0931 0.234
## traittarsus_s:genotypewt 0.240 0.3674 0.4400 0.5148 0.651
## traitSCT_s:genotypewt -0.274 -0.1543 -0.0945 -0.0294 0.114
## traitfemur_s:temp30 -0.933 -0.8698 -0.8389 -0.8025 -0.734
## traittibia_s:temp30 -0.919 -0.8505 -0.8076 -0.7698 -0.692
## traittarsus_s:temp30 -1.334 -1.2706 -1.2338 -1.2061 -1.147
## traitSCT_s:temp30 -0.958 -0.8774 -0.8368 -0.7999 -0.724
## traitfemur_s:genotypewt:temp30 0.220 0.3185 0.3608 0.4102 0.505
## traittibia_s:genotypewt:temp30 0.319 0.4229 0.4745 0.5277 0.626
## traittarsus_s:genotypewt:temp30 0.387 0.4720 0.5158 0.5593 0.648
## traitSCT_s:genotypewt:temp30 0.580 0.6750 0.7280 0.7773 0.877
And for the random effects.
summary(MMLM1.fit$VCV)
##
## Iterations = 2001:9991
## Thinning interval = 10
## Number of chains = 1
## Sample size per chain = 800
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## traitfemur_s:traitfemur_s.line 0.2965 0.0997 0.003524 0.003524
## traittibia_s:traitfemur_s.line 0.2494 0.0892 0.003155 0.003155
## traittarsus_s:traitfemur_s.line 0.2268 0.0842 0.002977 0.002977
## traitSCT_s:traitfemur_s.line 0.0297 0.0543 0.001921 0.001921
## traitfemur_s:traittibia_s.line 0.2494 0.0892 0.003155 0.003155
## traittibia_s:traittibia_s.line 0.2747 0.0911 0.003220 0.003220
## traittarsus_s:traittibia_s.line 0.2235 0.0805 0.002847 0.002847
## traitSCT_s:traittibia_s.line 0.0468 0.0523 0.001847 0.001847
## traitfemur_s:traittarsus_s.line 0.2268 0.0842 0.002977 0.002977
## traittibia_s:traittarsus_s.line 0.2235 0.0805 0.002847 0.002847
## traittarsus_s:traittarsus_s.line 0.2651 0.0880 0.003113 0.003113
## traitSCT_s:traittarsus_s.line 0.0805 0.0552 0.001951 0.001951
## traitfemur_s:traitSCT_s.line 0.0297 0.0543 0.001921 0.001921
## traittibia_s:traitSCT_s.line 0.0468 0.0523 0.001847 0.001847
## traittarsus_s:traitSCT_s.line 0.0805 0.0552 0.001951 0.001951
## traitSCT_s:traitSCT_s.line 0.1919 0.0621 0.002197 0.002568
## traitfemur_s:traitfemur_s.units 0.6441 0.0197 0.000698 0.000698
## traittibia_s:traitfemur_s.units 0.4014 0.0179 0.000634 0.000634
## traittarsus_s:traitfemur_s.units 0.2142 0.0150 0.000531 0.000453
## traitSCT_s:traitfemur_s.units 0.0257 0.0170 0.000599 0.000599
## traitfemur_s:traittibia_s.units 0.4014 0.0179 0.000634 0.000634
## traittibia_s:traittibia_s.units 0.6653 0.0219 0.000776 0.000776
## traittarsus_s:traittibia_s.units 0.1799 0.0144 0.000510 0.000482
## traitSCT_s:traittibia_s.units 0.0600 0.0173 0.000612 0.000576
## traitfemur_s:traittarsus_s.units 0.2142 0.0150 0.000531 0.000453
## traittibia_s:traittarsus_s.units 0.1799 0.0144 0.000510 0.000482
## traittarsus_s:traittarsus_s.units 0.5353 0.0177 0.000627 0.000627
## traitSCT_s:traittarsus_s.units 0.0847 0.0157 0.000555 0.000655
## traitfemur_s:traitSCT_s.units 0.0257 0.0170 0.000599 0.000599
## traittibia_s:traitSCT_s.units 0.0600 0.0173 0.000612 0.000576
## traittarsus_s:traitSCT_s.units 0.0847 0.0157 0.000555 0.000655
## traitSCT_s:traitSCT_s.units 0.7482 0.0243 0.000860 0.000860
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## traitfemur_s:traitfemur_s.line 0.16479 0.22595 0.2753 0.3457 0.5605
## traittibia_s:traitfemur_s.line 0.12555 0.18631 0.2322 0.2950 0.4862
## traittarsus_s:traitfemur_s.line 0.10937 0.16773 0.2065 0.2713 0.4333
## traitSCT_s:traitfemur_s.line -0.07347 -0.00275 0.0267 0.0617 0.1437
## traitfemur_s:traittibia_s.line 0.12555 0.18631 0.2322 0.2950 0.4862
## traittibia_s:traittibia_s.line 0.14585 0.20909 0.2569 0.3236 0.5014
## traittarsus_s:traittibia_s.line 0.11176 0.16435 0.2095 0.2645 0.4366
## traitSCT_s:traittibia_s.line -0.04938 0.01443 0.0424 0.0759 0.1565
## traitfemur_s:traittarsus_s.line 0.10937 0.16773 0.2065 0.2713 0.4333
## traittibia_s:traittarsus_s.line 0.11176 0.16435 0.2095 0.2645 0.4366
## traittarsus_s:traittarsus_s.line 0.14470 0.20043 0.2497 0.3094 0.4836
## traitSCT_s:traittarsus_s.line -0.00981 0.04533 0.0737 0.1091 0.2011
## traitfemur_s:traitSCT_s.line -0.07347 -0.00275 0.0267 0.0617 0.1437
## traittibia_s:traitSCT_s.line -0.04938 0.01443 0.0424 0.0759 0.1565
## traittarsus_s:traitSCT_s.line -0.00981 0.04533 0.0737 0.1091 0.2011
## traitSCT_s:traitSCT_s.line 0.10503 0.14808 0.1807 0.2246 0.3471
## traitfemur_s:traitfemur_s.units 0.60794 0.63077 0.6440 0.6567 0.6832
## traittibia_s:traitfemur_s.units 0.36780 0.38867 0.4016 0.4132 0.4387
## traittarsus_s:traitfemur_s.units 0.18585 0.20357 0.2142 0.2239 0.2435
## traitSCT_s:traitfemur_s.units -0.00857 0.01384 0.0256 0.0369 0.0601
## traitfemur_s:traittibia_s.units 0.36780 0.38867 0.4016 0.4132 0.4387
## traittibia_s:traittibia_s.units 0.62391 0.65005 0.6650 0.6792 0.7087
## traittarsus_s:traittibia_s.units 0.15161 0.17024 0.1799 0.1900 0.2095
## traitSCT_s:traittibia_s.units 0.02593 0.04848 0.0602 0.0715 0.0945
## traitfemur_s:traittarsus_s.units 0.18585 0.20357 0.2142 0.2239 0.2435
## traittibia_s:traittarsus_s.units 0.15161 0.17024 0.1799 0.1900 0.2095
## traittarsus_s:traittarsus_s.units 0.50501 0.52224 0.5345 0.5466 0.5723
## traitSCT_s:traittarsus_s.units 0.05495 0.07422 0.0847 0.0943 0.1188
## traitfemur_s:traitSCT_s.units -0.00857 0.01384 0.0256 0.0369 0.0601
## traittibia_s:traitSCT_s.units 0.02593 0.04848 0.0602 0.0715 0.0945
## traittarsus_s:traitSCT_s.units 0.05495 0.07422 0.0847 0.0943 0.1188
## traitSCT_s:traitSCT_s.units 0.70332 0.73160 0.7481 0.7636 0.7951
This is not the most friendly output, and it takes a while to get used to. However, we can see that we still have evidence for something interesting going, and we could extract the vectors of effects as we did above.
let’s look at the VCV matrix for the random effects of line in a slightly clearer way (as a matrix)
##' extract variance-covariance matrices for MCMCglmm objects
##' does not (yet) set cor, stddev, residual-var attributes
##' may be fragile: depends on group vars not having dots in them?
VarCorr.MCMCglmm <- function(object, ...) {
s <- summary(object$VCV)$statistics[,"Mean"]
grps <- gsub("^[^.]+\\.([[:alnum:]]+)$","\\1",names(s))
ss <- split(s,grps)
getVC <- function(x) {
nms <- gsub("^([^.]+)\\.[[:alnum:]]+$","\\1",names(x))
n <- length(nms)
L <- round(sqrt(n))
dimnms <- gsub("^([^:]+):.*$","\\1",nms[1:L])
return(matrix(x,dimnames=list(dimnms,dimnms),
nrow=L))
}
r <- setNames(lapply(ss,getVC),unique(grps))
return(r)
}
vv <- VarCorr(MMLM1.fit)
par(mfrow=c(1,2))
corrplot.mixed(cov2cor(vv$line),upper="ellipse")
corrplot.mixed(cov2cor(vv$units),upper="ellipse")
lme4 estimates (variables are in a different order)This is only a taste of what to do, and after this we could start asking questions about this genetic variance co-variance matrix.
Compare fixed effects:
tt <- tidy(MMLM1.fit)
tt2 <- tidy(lmer1)
tt_comb <- bind_rows(lmer=tt,MCMCglmm=tt2,.id="model") %>%
filter(effect=="fixed")
dwplot(tt_comb)+geom_vline(xintercept=0,lty=2)
glmmTMB, brms …FIXME:
allFit not working?)